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Abstract: This paper is concerned with global-in-time, nonoverlapping domain decomposition 
methods for the mixed formulation of the diffusion problem. Two approaches are considered: 
one uses the time-dependent Steklov-Poincare operator and the other uses Optimized Schwarz 
Waveform Relaxation (OSWR) based on Robin transmission conditions. For each method, a 
mixed formulation of an interface problem on the space-time interfaces between subdomains is 
derived, and different time grids are employed to adapt to different time scales in the subdomains. 
Demonstrations of the well-posedness of the subdomain problems involved in each method and a 
convergence proof of the OSWR algorithm are given for the mixed formulation. Numerical results 
for 2D problems with strong heterogeneities are presented to illustrate the performance of the two 
methods. 
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Methode de Decomposition de Domaine 
Espace- Temps pour les Problemes de Diffusion en 
Formulation Mixte 

Resume : Ce papier traite de methodes de decomposition de domaine sans re- 
couvrement, globales en temps, appliquees a un probleme de diffusion en formulation 
mixte. Deux approches sont considerees: l'une basee sur l'operateur de Steklov- 
Poincare, et l'autre basee sur une methode de relaxation d'onde optimisee (OSWR), 
utilisant des conditions de transmission de type Robin. Pour chaque methode, un 
probleme d'interface est ecrit sous forme mixte, les inconnues etant sur les interfaces 
espace-temps entre les sous-domaines, et differentes grilles en temps sont utilisees, 
adaptees aux differentes echelles temporelles dans les sous-domaines. Des demonstra- 
tions a la fois du caractere bien pose des problemes locaux dans les sous-domaines 
(intervenant dans chaque methode) et de la convergence de l'algorithm OSWR sont 
donnees en formulation mixte. Des resultats numeriques pour des problemes 2D avec 
de fortes heterogeneites sont presentes pour illustrer les performances des deux meth- 
odes. 

Mots-cles : formulation mixte, decomposition de domaine espace-temps, prob- 
leme de diffusion, operateur de Steklov-Poincare dependant du temps, methode de 
relaxation d'onde optimisee, maillages non-conformes en temps 
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1. Introduction. In many simulations of time-dependent physical phenomena, 
the domain of calculation is actually a union of several subdomains with different 
physical properties and in which the time scales may be very different. In partic- 
ular, this is the case for the simulation of contaminant transport around a nuclear 
waste repository, where the time scales vary over several orders of magnitude due to 
changes in the hydrogeological properties of the various geological layers involved in 
the simulation. Consequently, it is inefficient to use a single time step throughout 
the entire domain. The aim of this article is to investigate, in the context of mixed 
finite elements [5j I30| , two global-in-time domain decomposition methods well-suited 
to nonmatching time grids. Advantages of mixed methods include their mass conser- 
vation property and a natural way to handle heterogeneous and anisotropic diffusion 
tensors. 

The first method is a global-in-time substructuring method which uses a Steklov- 
Poincare type operator. For stationary problems, this kind of method (see [29i [341 [28] ) 
is known to be efficient for problems with strong heterogeneity. It uses the so-called 
Balancing Domain Decomposition (BDD) preconditioner introduced and analyzed in 
[2"4l, 125] , and in [5] for mixed finite elements. In brief, the method "involves at each 
iteration the solution of a local problem with Dirichlet data, a local problem with 
Neumann data and a "coarse grid" problem to propagate information globally and to 
insure the consistency of the Neumann problems" [5]. 

The second method uses the Optimized Schwarz Waveform Relaxation (OSWR) 
approach. The OSWR algorithm is an iterative method that computes in the sub- 
domains over the whole time interval, exchanging space-time boundary data through 
more general (Robin or Ventcel) transmission operators in which coefficients can be 
optimized to improve convergence rates. Introduced for parabolic and hyperbolic 
problems in [9], it was extended to advection- reaction-diffusion problems with con- 
stant coefficients in [26]. The optimization of the Robin (or Ventcel) parameters was 
analyzed in |1Q[ [T] and extended to discontinuous coefficients in [TTl 12]. Extensions 
to heterogeneous problems and non-matching time grids were introduced in [111 15]. 
More precisely, in [3j[T5], discontinuous Galerkin (DG) for the time discretization of 
the OSWR was introduced to handle non-conforming time grids, in one dimension 
with discontinuous coefficients. This approach was extended to the bidimensional 
case in [171 02] • One of the advantages of the DG method in time is that a rigorous 
analysis can be carried out for any degree of accuracy and local time-stepping, with 
different time steps in different subdomains (see |17l I18| ). A suitable time projection 
between subdomains was obtained by an optimal projection algorithm without any 
additional grid, as in [T5J. These papers use Lagrange finite elements. An extension 
to vertex-centered finite volume schemes and nonlinear problems is given in |14| . The 
classical Schwarz algorithm for stationary problems with mixed finite elements was 
analyzed in [7J. 

In this work, we extend the first method to the case of unsteady problems and 
construct the time-dependent Steklov-Poincare operator. For parabolic problems, we 
need only the Neumann-Neumann preconditioner |22] as there are no difficulties con- 
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cerning consistency for time-dependent Neumann problems. Of course one could make 
use of the idea of the "coarse grid" to ensure a convergence rate independent of the 
number of subdomains. However, we haven't developed this idea here. The conver- 
gence of a Jacobi iteration for the primal formulation is independently introduced and 
analyzed in |21| . 

For the second method, an extension of the OSWR method with Robin trans- 
mission conditions to the mixed formulation is studied and a proof of convergence is 
given. For each method a mixed formulation of an interface problem on the space- 
time interfaces between subdomains is derived. The well-posedness of the subdomain 
problems involved in the first approach is addressed in [531 EL through a Galerkin 
method and suitable a priori estimates. In this paper we present a more detailed ver- 
sion of the proof for Dirichlet and extend these results to to prove the well-posedness 
of the Robin subdomain problems involved in the OSWR approach. In [311152] demon- 
strations using semigroups are given for nonlinear evolution problems. For strongly 
heterogeneous problems, it is natural to use different time steps in different subdo- 
mains and we apply the projection algorithm in [13] adapted to time discretizations 
to exchange information on the space-time interfaces, for the lowest order DG method 
in time. We show the numerical behaviour of both methods for different test cases 
suggested by ANDRA for the simulation of underground nuclear waste storage. A 
preliminary version of this work was given in |19| . 

The remainder of this paper is organized as follows: in the next section we present 
the model problem in a mixed formulation. We prove its well-posedness for Dirichlet 
and Robin boundary conditions in Section[3J In Sectional we introduce the equivalent 
multidomain problem using nonoverlapping domain decomposition and describe the 
two solution methods. A convergence proof for the OSWR algorithm for the mixed 
formulation is given. In Section [5l we consider the semi-discrete problems in time 
using different time grids in the subdomains. In section [6j results of 2D numerical 
experiments showing that the methods preserve the order of the global scheme are 
discussed. 

2. A model problem. In this section we define our model problem and show 
the existence and uniqueness of its solution. For an open, bounded domain fi of 
R d (d = 2, 3) with Lipschitz boundary dfl and some fixed time T > 0, we consider 
the following time-dependent diffusion problem 

wd t c+ V- (-Z>Vc) = /, inOx(0,T), (2.1) 

with boundary and initial conditions 

c = 0, oncK7 x (0,T), 
c(-,0)=c o , inO. (2.2) 

Here c is the concentration of a contaminant dissolved in a fluid, / the source term, 
lo the porosity and D a symmetric time independent diffusion tensor. We assume 
that uj is bounded above and below by positive constants, < u;_ < lj{x) < w+, 
and that there exists 5- and 5+ positive constants such that £ T ,D _1 (x)£ > £_|£| 2 , 
and |D(x)^| < for a.e. x € Q and V£ € K d . For simplicity, we have imposed a 

homogeneous Dirichlet boundary condition on d£l. In practice, we may use non- 
homogeneous Dirichlet and Neumann boundary conditions for which the analysis 
remains valid (see Section [3] for the extension to Robin boundary conditions). 
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We now rewrite (|2.ip in an equivalent mixed form by introducing the vector field 
r := -DVc. This yields 



u;<9 t c + V-r =/, infix(0,T) 
Vc + XrV =0, inOx(0,T) 



-1_ n ;„o.Wn^ (2-3) 



To write the variational formulation for (|2.3p (see 0[3O]), we introduce the spaces 

M = L 2 (fi) and S = H (div, ft) . 

We multiply the first and second equations in (|2.3|) by /x € A'/ and w € £ respectively, 
then integrate over ft and apply Green's formula to obtain: 

For a.e. t € (0,T), find c(i) € M and r(t) G £ such that 

— (wc,/i) + (V -r,/i) =(f,fx), V/USM, / 2 4 x 

-(V-i>,c) + (D~V,w) =0, Vv G S, 



together with initial condition (12.21) . 

Here and in the following, we will use the convention that if V is a space of 
functions, then we write V for a space of vector functions having each component 
in V. We also denote by (-,-) the inner product in L 2 (Q) or Z 2 (ft) and || • || the 
L 2 (ft)-norm or L 2 (f2)-norm. 

The well-posedness of problem (|2.4[) is shown in in [2"31 H], with an argument based 
on a Galerkin's method and a priori estimates: 

Theorem 2.1. If f is in L 2 (0, T; £ 2 (ft)) andc mi^ft) then problem ([231, 
has a unique solution 

(c,r) £H\0,T;L 2 (n)) x (i 2 (0, T; F(div, 0)) n i°°(0, T; L 2 (n))) . 

Moreover, i/D is m l^ 1,00 (n), / in H 1 ^, T; L 2 (fl)) and c in H 2 (fl) n Fo(^) </lerl 

(c,r) Gl^ 1 ' oo (0,T;L 2 (f2)) x (i°°(0, T; £T(div, ft)) n ff 1 (0,T;L 2 (f2))) . 

Remark. We give the proof of Theorem thrml in the finite dimensional setting 
since some technical points (those involving dtr, or r at time t = 0) can only be 
defined by their finite dimensional Galerkin approximation. This is not surprising 
given the differential-algebraic structure of system (|3.2[) : the second equation has no 
time derivative. In DAE theory it is well known that the algebraic equations have 
to be differentiated a number of times (this is what defines the index) , and that this 
imposes compatibility conditions between the initial data (note that r(0) is not given). 
The index has been extended to PDEs, see for instance |27j . 

The proof of Theorem l2.1l is carried out in several steps: in Lemma [2~2l we first con- 
struct solutions of certain finite-dimensional approximations of (|2.4p , then we derive 
suitable energy estimates in Lemma 12.31 and prove the first part of the theorem. The 
higher regularity of the solution is obtained from the estimates given in Lemma 12.41 

We need first to introduce some notations: Let {/i„ | n G N} be a Hilbert basis 
of M and {v n | n G N} be a Hilbert basis of E. For each pair of positive integers n 
and m, we denote by M n the finite dimensional subspace spanned by {/ii}" =1 , and 
S m the finite dimensional subspace spanned by {Vi}^l 1 . Now let c„ : [0,T] — > M n 
and r m : [0, T] — > E m be the solution of the following problem 

(wdtCn, fM) + (V • r m ,fii) =(f(t), Hi), Vi = l,...,n, , . 

-(V-t;„c„) + (D- 1 r m ,t;,) =0, Vj = 1,...,to, [ -> 



(i 



with 

(c„,(0),Hi) = (c ,/ii), Vi = l,...,n. (2.6) 

Lemma 2.2. (Construction of approximate solutions) For each pair (n, m) £ N 2 , 
n,m > 1, i/iere exists a unique solution (c n ,r m ) to problem (]2 . 5[) . 
Proof. We introduce the following notations 

(F n (t))i = (f(t),fj,i), (C )» = (co,j"t), (W" n )ij = (b}Hj,m), VI < i,j < n, 

(A m )y = (D~ 1 v j ,v i ) 7 VI < i,j < m, (B nm )ij = (V -v^m), V 1 < i < n, 1 < j < m. 

We also denote by C„ (t) the vector of degrees of freedom of c„ (i) with respect to the 
basis and R m (t) that of r m (t) with respect to the basis {wi}7=i- With this 

notation, (12. 5p may be rewritten as 

W„^(t) + B nm R m (t) = F„(i), (2.7a) 

-B T nm C n {t) + A m R m {t) = 0, (2.7b) 
C n (0)=C o . (2.7c) 

As A m is a symmetric and positive definite square matrix of size m (because of the 
assumptions concerning D), A m is invertible. Thus (|2.7b[) implies 

R m {t) = A^Bl m C n {t). (2.8) 

Substituting (|2.8j) into ()2.7a[) and as W"„ is invertible, we obtain 

dC n , 



(It 



■(*) + W/^jjBLCl'l = ^(t), for a.e. t e [0,T]. (2.9) 



This is a system of n linear ODEs of order 1 with initial condition (I2.7cp . Hence, 
there exists a unique function C„ € (C([0,T])) n with —y— € (L 2 (0 1 T)) n satisfy- 
ing (HU) and (j2~Tcj) (see 0). From we obtain R m £ (C([0,T])) m such that 
dR 

m G (L 2 (0,T)) and then (c n ,r m ), which is the unique solution to (|2.5[) . □ 

In the next step, we derive some suitable a priori estimates similar to those given 
in [23J but in a more detailed manner. 

Lemma 2.3. There exists a constant C independent of n and m such that 

ll c n||,L°°(0 : T;L 2 (Sl)) + \\dtC n \\ L 2 (0,T;L 2 (n)) + \\ r m \\ (0,T;L 2 (f2)) + ||»"m||i»(0,T;H(div,n)) 

< C(ll c o||ffi(o) + II/IIl2(o,T;L2(o))), Vn,m > 1. 

Proof. We prove this lemma by deriving successively the estimates on c n , dtc n 
and r m , and finally on V • r m for the i?(div, f2)-norm. 

• Let n, m > 1 and take c„(t) e M n and r m (t) 6 E m as the test functions in (|2.5[) 

(ujd t c n , c n ) + (V • r m , C n ) =(/,c„), 
-(V ■ r m ,Cn) + (D _1 r m ,r m ) =0. 

Adding these two equations, we obtain 

(ujd t c n ,c n ) + (D~ 1 r m ,r rn ) = (/,c„). 
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Using the properties of lo and D, and applying the Cauchy-Schwarz inequality, we get 

(ujd t c n ,c n ) = ^^-(uc n (t),c n (t)) > 'Y^\\ c "( t )W 2 ^ 
(ZrV m (t),r m (i)) > S-\\r m {t)\\, 

(f(t),C n (t)) < \\f(t)\\ \\ Cn (t)\\ < 2-^ll/WH 2 + ^l|Cn(i)|| 2 . 

As w_ > 0, we deduce that 



^llcnWII 2 + —\\r m (t)\\ 2 < -hr\\f(t)\\ 2 + \K(t)\\ 2 



Integrating this inequality over (0,i) for t £ [0,T], we find 

Ir / 

U! 



\K(t)\\ 2 + ^ J\\r m (s)\\ 2 ds < ||c(0)|| 2 + -i- / /(.s)||-</.s+ / r„ { s)f<l.s. (2.1(1) 



J - Jo w - Jo Jo 



since ||c„(0)|| 2 = ^(c ,/^) 2 < ^(c ,^) 2 - ||cq| 



^0, fMJ — ||^0|' 2 
i=l i=l 

Thus (I2J0)) implies 

||C„W|| 2 < (||C0|| 2 + ^||/||! 2( 0,T ; ^(U))) + [ \\Cn(s)\\ 2 d S . 







Applying Gronwall's lemma, there exists C independent of n or m such that 

<C(|| C o|| 2 + ||/|| 2 L2(0 , T;i2(f2 ))), (2-11) 

• Now we derive the estimate for dtc n : Taking 3 t c„ £ M n as the test function in the 
first equation of (|2.5p . we obtain 

(ujd t c n ,d t c n ) + (V ■ r m ,d t c n ) = {f,d t c n ). (2.12) 

Differentiating the second equation of (|2.5[) with respect to t, we find 

-(V ■v,d t c n ) + (D- 1 d t r m ,v)^0, VueS m . (2.13) 

Then we take r m as the test function in (|2.13l) 

{D- l d t r mi r m ) - (V T m ,a t c„) = 0. (2.14) 

Adding (|2~T2]) and ([2TT4]) . we see that 

(ujdtc n ,d t c n ) + (D~ 1 d t r m ,r m ) = (f,d t c n ). 

As D is symmetric and positive definite, by applying the Cauchy-Schwarz inequality 
to the right hand side as well as using the property of w, we obtain 

uj_\\d t c n (t)\\ 2 + ^-\\VD^r m (t)\\ 2 < — ll/WII 2 - (2-15) 

at uj- 
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Integrating (|2~T5j) over (0,t) for t € [0,T], we find 

u;_ / ||cW(s)|| 2 ds+||V^r m (f)|| 2 < ||VS rT r m (0)|| 2 + — / 
Jo w - Jo 



||/(a)|| 2 ds. (2.16) 



To bound ||r m (0)|j, we take r m 6 S m as the test function in the second equation 
of and let t = 

(irV^r^O)) = (V -r m (0),c„(0)). (2.17) 

Noting that (|2.17l) holds for all n,m > 1, we bound the left-hand side as before and 
let n — > oo. Since c ra (0) — > Cq in L 2 (f2) and Cq G iJg(^)> we have by Green's formula 

<5-IMO)|| 2 < (V-r ro (0),c ) = (r m (0),-Vc ) < ||r w (0)|| ||Vc ||. 

Thus 

||r m (0)|| < C\\co\\ Hi[n y (2.18) 

This along with (|2~16|) yields 

l|c , t c n|lL 2 (0.T;L 2 (O)) + ll r "i|li°°(0.T;L 2 (f2)) ^ C(ll c o|lfl-i(f2) + II / 1 1 Z,2( ,T;i 2 (fl) ) ) ' VVl, m > 1. 

(2-19) 

There only remains to show that ||V • fr"m||.L 2 (o,T;i, 2 (fi)) i s bounded. 
• Fixing m > 1, as V • r m (t) £ M we can write 

oo 

V -r m (i) = for a - e - f e (0,T), (2.20) 

1=1 

where = (V ■ r m (t),/ij). Now we fix n > 1 and multiply the first equation 

of (|2.5p by sum over « = 1, we see that 



(VT^f^w) < i(||/|| + q|a tC „||) 2 + i|| X)4wf. (2-21) 

i=l i=l 

Integrating with respect to time and recalling (|2.19p , we find 
/ (V-r ro ,^e>i)^<C(||/||l 2(0 , X;i2(n)) + 110011^^) + - / || 2 toll 2 *- 

Jo i=l Jo i=l 



Let n — > oo and recall (|2.20p . we obtain 

||V • r m \\ 2 dt < C(||/|| 2 2(0 T . i2(n)) + || co || 2 Hi{Q) ) + i / ||V • r m || 2 rft. 
Thus 

ll^ 7 ' r ™lli 2 (0,T;L 2 (n)) ^ C(ll/lli 2 (0,T;L 2 (n)) + 1 1 c II iyi (fi) ) ■ 

On the other hand, by recalling inequality (|2.10p with t — T and by (|2.1ip . we find 

ll r "i|l,L 2 (0,T;L 2 (n)) 

<C'(|| Co || 2 + ii/n 

L 2 (0,T;L 2 (n)))- 
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Hence, 



l r '"llL 2 (0,T;H(div,fi)) — ll r "i|ll 2 (0,T;(L 2 (n)) 2 ) + ll^ 7 ' r "i|ll 2 (0,T;L 2 (n)) 



< C(I!/IIl 2 (0,T;L 2 (O)) + ll c 0||ffi(O))> Vm > 1 



which ends the proof of Lemma 12.31 □ 

We now prove the first part of Theorem 12.11 there exists a unique solution (c, r) 
in iT^O.TjL 2 ^)) x L 2 (0, T: H(div, f2)) n L°°(0, T; L 2 (f2)) of problem (|2~5)) . 

Proof. The proof of the first part of Theorem 12.11 follows the following steps:. 

• Lemma l2~3"l implies that for the sequences {c n }^ =1 and {r m }^ =1 defined by (|2.5|) and 
(EH), {cn)n=i is bounded in L 2 (0, T; L 2 (Q)), {d t c n }™ =1 is bounded in L 2 (0, T; L 2 (VL)) 
and {r m }^ =1 is bounded in L 2 (0, T; ff(div, Q)) n L°°(0, T; L 2 (0)). Thus, there exist 
subsequences, still denoted by {cn}^^ and {r m }m=i an< ^ functions c € L 2 (0, T; £ 2 (£!)) 
with <9 t c G L 2 (0,T ;J L 2 (ft)) and r G L 2 (0, T; iJ(div, O)) n L°°(0, T; L 2 (f2)) such that 

c„ -cin L 2 (0,T ;J L 2 (fi)), 

a t c„ -5 t cinL 2 (0,T;L 2 (r!)), (2.22) 
r m ^rin L 2 (0, T; i?(div, fi)). 

• Next let jj G ^([0, T]; M na ), w G ^([0, T]; E mo ) for n ,m > 1. We choose n > n 
and m > mo, take 77 and iy as the test functions in (|2.5I) and then integrate with 
respect to time 



T /-T 

(ujd t c n , 77) + (V • r m ,r])dt = / (/, 77)*, 



/ -(V -w,cn) + (D- 1 r m ,w)dt =0. 
Jo 

Because of the weak convergence in ()2.22j) . we also have 



(2.23) 



M t c,77) + (V-r,77)^ = / (/, 77)dt, 

/ -(V ■w,c) + (D- 1 r,w)dt =0. 
Jo 



(2.24) 



Since the spaces of test functions r],w are dense in L 2 (0,T; M) and L 2 (0,T;E) re- 
spectively, it follows from |2~24| that (f2T4|) holds for a.e. t G (0,T) (see (8]). 

• There remains to show that c(0) = cq. Toward this end, we take 77 G C ([0, T]; M no ) 
with ??(T) = 0. It follows from the first equation of (|2.24p that 



(ud t r,,c) + (V-r,r,)dt= / (/,r?)di + (wc(0),»j(0)). (2.25) 

Jo 

Similarly, from the first equation of (|2.23[1 we deduce 

,T 

(ud t r),c n ) + (V -r m ,r))dt = / (/, 77)* + (wc„(0), 77(0)). 



Using (|2.22p . we obtain 



(ijd t r],c) + (V-r,r))dt= f (/, 77)^ + (cjc , 77(0)), (2.26) 
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since c„(0) — > Co in L 2 {VL). As 77(0) is arbitrary, by comparing (|2.25p and (|2.26[) we 
conclude that c(0) = cq. 

• For the uniqueness, as the equations are linear, it suffices to check that c = and 
r = for / = and cq = 0. To prove this, we set fi = c and v = r in \2A\ (for / = 0) 
and add the two resulting equations: 

~(ojc,c) + (D- l r,r)=0. 

Using the property of a; and the fact that (D~ 1 r : r) > S-\\r\\ 2 > 0, then integrating 
with respect to t we see that 

uj_\\c(t)\\ 2 + 25- / \\r(s)\\h m d S < 0, fora.e. t£ (0,T), 
Jo 

where c(0) = c = 0. Thus c = and r = for a.e. t € (0,T). □ 

We now prove the second part of Theorem 12.11 The higher regularity of the 

solution to (|2.3p is obtained by using the following lemma. 

Lemma 2.4. (Estimates for improved regularity) Assume that D is in W 1 '°°{Q), 

c mH 2 (n) n #o(n) and fin H 1 {Q,T\L 2 {Q)) then 

\\dtc\\L^(0,T-L 2 (n)) + ll r lli°=(0,T;ff(div,a)) + ||%"||£2(0,T;L 2 (n)) 
< C(ll/l|H 1 (0,T;L 2 (f2)) + ||co||H 2 (n))- 

Proof. As / e J ff 1 (0,T; J L 2 (f7)), the solutions of the ODE system flU} are more 
regular in time than before (i.e. up to second-order time derivatives). 

Let n,m> 1. First, we differentiate the first equation of (|2.5p with respect to t 

(udttCn, m) + (V • dtr m ,(H) = (9t/,M»)j Vi = 1, ...,n, 

then we take 9tc n as the test function 

(ujdttc n ,dtc n ) + (V • d t r m ,d t c n ) = (d t f, d t c n ). (2.27) 

Similarly, we differentiate the second equation of (|2.5[) with respect to t 

(D- l dtr m ,Vi) - (V • Vi,d t c n ) = 0, Vi = l,...,m, 

and take <9tr m as the test function 

(D _1 5 t r m , %Vi) ~ (V • d t r m ,d tCn ) = 0. (2.28) 

Adding (f2~27j) and (|2~28|) . we find 

(ujdttc„,d t c n ) + (D~ 1 d t r m ,d t r m ) = {dtf, d t c n ). 

Bounding (D~ l d t r m , d t r m ) > <5_ ||<9 t r m || 2 , using the assumption about uj and applying 
the Cauchy-Schwarz inequality, we obtain 

^\\d t c n \\ 2 + — \\d t r m \\ 2 < -^\\d t f\\ 2 + \\d t c n \\ 2 . 
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For each t £ [0,T], we may integrate over (0,t) to obtain 

\\d t c n (t)\\ 2 + — I \\dtr m \\ 2 ds<\\d t c n (0)\\ 2 + \ ( \\d t f\\ 2 ds+f \\d t c n \\ 2 dd,2.2S)) 
Jo u_ j a j Q 

In order to bound ||c?tc n (0)||, we use the first equation of (|2.5[) (with 3 t c„ as the test 
function, at t — 0) to obtain 

||9 tCn (0)||<C(||VT ro (0)|| + ||/(0)||). 

Using the second equation of ()2.5p at t — 0, and then let n — > oo to get 

(£rV m (0) + Vco,v) =0, V» G E m . 

Thus, using density argument and cq £ H^(n)nH 2 (fl), we obtain D~ 1 r m (0) = — Vco 
in iJ x (fi). Then, we bound 

||a 4Cll (0)|| 2 <C(|| C0 ||^ (o) + ||/(0)|| 2 ). (2.30) 

Replacing ([2301) in (f2~2"§)) , we obtain 

\\d t c n (t)\\ 2 + — f \\d t r m fds 

<C(|| C0 || 2 ff2(o) + ||/|| 2 ffl(0iT:L2(o)) )+ / ||5 t c„|| 2 d S . (2.31) 

Jo 

It now follows from ()2.31l) and Gronwall's lemma that 
2S 

\\dtC n \\ L°-=(0,T;L2(Q)) + ~ H^t r m|lL 2 (0,T;L 2 (f2)) < C(ll c II ff2(n) + || /|| H i (0,T:L 2 (f2))X 2 -32) 

Recalling (I2.2ip and using (|2.32p . we obtain 

(VT m ,^Cft) < C(\\c \\ 2 HHn) + \\f\\ 2 H i {0 . T . L 2 m ) + gUj^&wH 2 . 
t=l i=l 

Then, let n — > oo, we see that 

II V • r m \\ 2 LO a( T . L 2( n ^ < C(\\c \\ 2 H 2^ + ll/llfl-i(o,T ; L 2 (n)))- 
This along with (|2.19p gives 

||*"m||ioo(o,T;H(div,f2)) - ^(11 c II H 2 (n) + II / II H 1 (0,T:L 2 (J2)) ) ■ (2.33) 

The lemma now follows from (pT32|) . (f2~33]) and (|2~22|) . □ 

In the sequel, we will consider two domain decomposition methods for solving 
(|2.4p . ()2.2p . The first one involves local Dirichlet subproblems whose well-posedness 
is an extension of Theorem 12.11 In the second approach, the optimized Schwarz 
waveform relaxation method, we shall impose Robin transmission conditions on the 
interfaces. Thus, we extend the well-posedness results above to the case of Robin 
boundary conditions. 
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3. A local problem with Robin boundary conditions. In this section, we 
consider problem (|2.1[) - (|2.2p with Robin boundary conditions on dft x (0,T) : 

-r-n + ac = g, on dft x (0, T), (3.1) 

where a defined on dQ is a time independent positive, bounded coefficient and g is a 

space-time function. We define a := — and suppose that < k± < a < Ki a.e. in dfl. 

a 

We denote by (•, and |j • \\qq the inner product and norm in L 2 (dVl) respectively. 
To derive a variational formulation corresponding to boundary condition (|3.ip . we 
introduce the following Hilbert space 

E = H(dW, Q) := {v e H(div, Q)\ v ■ n e L 2 {dn)}, 

equipped with the norm 

IHw(div,n) : = ll«lk(div,n) + ||» • »||eo- 
The weak problem with Robin boundary conditions may now be written as follows: 

For a.e. t € (0,T), find c{t) e M and r(t) € E such that 

(o;9 t c,M) + (VT,/x) =(/,M). V/xeM, 
-(V • v, c) + (-D~ 1 r,w) + (or • n,u • n)asi = -~{ag,v ■ n)dn, Vu e E. 

Theorem 3.1. /// is in L 2 (0, T; L 2 (fl)), g in ff^OjT; L 2 (dVL)) andc inH 1 ^), 
then problem (|3.2p . f|2 . 2|) /ias a unique solution 



(c,r)G H 1 {0,T;L 2 {n))x (L 2 (0, T; "H(div, 17)) n L°°(0, T; L 2 (f2))) . 
Moreover, j/D is in W ll00 (fl), / m H 1 [Q,T] L 2 {VL)) and c m iJ 2 (^) f/iera 

(c,r)e ^'"(O.T;!^!))) x (i oo (0,T;-H(div,fi))ni? 1 (0,r;i 2 (O))) . 

Proof. The proof of Theorem l3 . 1 1 relies on energy estimates and Gronwall's lemma, 
together with a Galerkin method, as for the proof of Theorem 12.11 We only present 
here parts of the proof that are different from those of the proof of Theorem 12.11 We 
construct the finite-dimensional approximation problems to (|3.2p as follows 



(ud t Cn, Mi) + (V • r m , fit) =(/,Mi)) 
-(V -Vj,c n ) + (D~ l r m ,ij) + (dr m n,«j -n),^ = (-<%,«.; -n^n,! < j < m, 



where c„ € M n , r m € E m and = 1, . . . , m is the basis of E m . We then rewrite 
(|3.3p in matrix form as in (|2.7p : 

dC 

W n —^(t) + B nm R m (t) = F n (t), 
at 

~BnmPn{t) + -^m-RmW = G m (t), 

where .R m is the vector of degrees of freedom of r m with respect to the basis {vi}™^^, 
(A m )ij = (D^v^Vi) + (avj n,Vi- n) d n, VI <i,j < m, 
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is symmetric and positive-definite, 

(B nm )ij = (V -Vj,fjbi) and (G m {t)) l = (-ag(t),Vi -n)gn, V 1 < i < n, 1 < j < m. 

Thus, there exists a unique solution (c n ,r m ) to (|3.3|l . 

Now to prove the existence of a solution to (|3.2|) . we derive suitable energy esti- 
mates in the same manner as in Section[2]but with an extra term r n on the boundary. 

Lemma 3.2. Let f e L 2 (0, T; L 2 (VL)) , g G i/ 1 (0, T; L 2 (dQ.) and c € ff 1 ^). 
TTie following estimates hold 

(*) II c IIl~(0,T;L 2 (S1)) + lk|| L 2( ,T;L 2 (n)) + ll r ■ n IU 2 (0,T;L2(en)) 

< C(II c o||l 2 (I2) + ||/||L 2 (0,T;L 2 (f2)) + \\g\\L 2 (0,T;L 2 (dn))), 

(it) \\d t c\\ L 2( T . L 2 {n)) + ||r|| iTC(0iT . L 2 (n)) + \\r ■n|| L °o( 0iT;L 2 (an) ) 

< CdlcolliJ^n) + ll/IU 2 (0,T;L 2 (O)) + ||5l|ff 1 (0,T;L 2 (aO))) 1 

(jn) ||r||ia(o,T;«(dOT,n)) < C(||co||ifi(n) + I!/IU 2 (0,T;L 2 (£1)) + ||.9||//i(0,T:L 2 (9Q)))- 

Lemma 3.3. (Estimates with greater regularity) Assume that D is in W 1,00 (S}), 
c m H 2 (n), f in ff 1 (0, T; £ 2 (S1)) and 3 m H x (0, T; L 2 (dfl)), then 

\\dtc\\L°°(0,T;L*(n)) + ll r llL~(0,T;H(div,O)) + ll^ r llL 2 (o,T;L 2 (n)) 
< ^(ll/lli? 1 (0,T;i 2 (Q)) + 1 1 c o||// 2 (S3) + ll5ll// 1 (0,T;L 2 (an)))- 

Proof, (of Lemma \3.2\l . In order to prove (i), as before, we take c„ and r m as 
test functions in (|3.3|) and add the two equations: 

(u>d t c n ,c n ) + (D~ 1 r m ,r m ) + (or m n,J- m ■ n) m = (f,c n ) + (-ag,r m ■ n) gQ . 

The assumptions concerning u>, D and a give 

(ud t c n ,c n ) > -y^ll c «ll 2 > (D^ 1 r rn ,r m ) > <S_||r m || 2 , (ar m n,r m - n) dn > Ki||r m -n||| n , 
and the Cauchy-Schwarz inequality: 

|(/, C n)|<||/||||c„||<^||/|| 2 + ^|| C „|| 2 . (3.4) 
Similarly, for each e > 0, 

| - (ag,r m -n) dn \< K 2 \\g\\dn ||r ro -n|| a n < n 2 (^\\g\\g Q + |||r TO • nf d ^j . (3.5) 
Choosing e = — , we then obtain 
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Integrating this inequality over (0, t) for t € (0, T], and using ||c„(0)|| < ||co|| , we get 

||c„(t)|| 2 + — [ t \\r m (a)\\*da + £- fwrnis). n\\% n ds 
w - Jo w - Jo 

< C (\\c \\ 2 + ||/|||2 (0iT;L 2 (o)) + \\g\\h(0,T;L2(dn))) + / \\Cn (s) f ds, 



1 K 



with C — max(l, —5-, ). Then an application of Gronwall's lemma completes 

UJ_ LO-Ki 

the proof of (i). 

For (ii), we follow the same steps as in (|2.12|) - (|2.15|) : taking d t c n € L 2 (0,T; M) 
as the test function in the first equation of (|3.3I) . we obtain 



(wftc, d t c) + (V ■ r m , d t c) = (/, d t c). (3.6) 

Differentiating the second equation of (|3.3|) with respect to t, we obtain 

-(V • v,d t c n ) + {D~ 1 d t r m ,v) + {ad t r m ■ n,v ■ n) dn = -(&d t g,v- n) 9 n, Vv € E. 

Then we take u = r m in the previous equation and add the resulting equation to (|3.6p 
to obtain 

{ujd t c n7 d t c n ) + {D~ l d t r m ,r m ) + {ad t r m -n,r rn -n) dn = {f,d t c n ) - (ad t g,r m ■ n) gn . 

As D is symmetric and positive definite, by applying the Cauchy-Schwarz inequality 
to the right hand side as well as using the property of w, we obtain 

w - ||5tc| ' 2 + \jt llV ^ r ' j2 + YJt ]]rm ' n|lln - 1 iLdtc) 1 + 1 Mt3 ' rm ' n)dn 1 ■ 

We then apply the Cauchy-Schwarz inequality for the right-hand side (as in (|3.4[) . (|3.5p , 
replacing c and g by <9 t c and 9j<?), and take e = — , C — max( — , — ) to obtain 

u-\\dtcn\\ 2 + j t \\ VD^r m \\ 2 + Kl |||r m -n||| n < C (||/|| 2 + \\d t g\\% a ) + Kl \\r m -nf dn . 
Integrating over (0, t) for t G [0, T], we find 



w- / ||9 t c„(s)|| 2 ds + ||\/^ T r m (<)|| 2 + Kl ||r m (t)-n||^ 
Jo 

< C(\\f\\h(p,T;W(p)) + \\9\\m(o.T-LHdn))) + II ^D ZT r m (0) || 2 + Kl \\r m (0) • n||| n 



+ «i / ||r TO (s) -nlUnds. (3.7) 
Jo 



So there only remains to bound the term (||v'£)~ 1 r m (0) || 2 + ki ||r m (0)-fl|||n)- Toward 
this end, we use the second equation of (I3.3|) with v — r m and for t = to obtain: 



<5_||r m (0)|| 2 + Ki||r m (0) -n|| 2 < (V • r m (0), c„(0)) + (-d ff (0),r m (0) • n) an . 
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Let n — > oo, as c„(0) — > cq we have 

<5_||r m (0)|| 2 + Ki||r m (0) -n\\ 2 < (V-r m (0),co) + (-d S (0),r m (0) ■ n) 9n 

< (-r m (0),Vc ) + (c - agr(0),r m (0) -n) an 

< y l|r ro (0)f + 2T ll Vc °H 2 + flM ) ' "II 2 + ^|H C " - 9(0)\\an, 



<5_||r m (0)|| 2 + Ki||r m (0) -n|| 2 < C (iN^n) + \\g\\ 

This along with (|3.7j) and Gronwall's lemma yields (it). We now estimate ||V • r m || 2 
as in Section [2] we derive (|2 . 2 1 1) from (|2.20p and the first equation of (|3.3j) (after 
multiplying by £' n (i) and summing over i = 1, ...,n). Then, using the bound for 

ll^c|| L 2( 0iT . i 2 (o)) in (it), we obtain 

< C(\\co\\ + 11/11 L 2 (0,T;L 2 (fi)) + IMIi?i(0,T;.L 2 (afi)))' (3-8) 

This along with (i) gives 

H r ™llL 2 (0,T;-H(div,f2)) 

<C(|| c o||_ffi(f2) + 11/11 L 2 (0,T;L 2 (n)) + IMI 



and the proof of Lemma 13.21 is completed. □ 

We now prove Lemma 13.31 for the higher regularity of the solution to (|2.3| . 

Proof. ( of Lemma Iff. 3)) . Let n, m > 1. Differentiate both equations of (|3.3| 
with respect to t, take // = 9tc„ and v = 9*r m as the test functions, and add the two 
resulting equations to obtain 

(ujdttc n ,dtc n ) + (D _1 5 t r m , d t r m ) + (ad t r m ■ n, d t r 

= (d t f, d t c n ) - (ad t g, d t r m ■ n) gn . 

Then, the assumptions concerning lj, D, a, and the Cauchy-Schwarz inequality give 

^^ ll ^ C - ||2+f5 - |l ^ T *^ l|2 + ^ L||atT "-- iri11 ^ ^ ^^H^/l^ + ^-ll^lll^ + ^ll^^ll 2 - 
Integrating this inequality over (0,i), for for t £ (0, T], we obtain 

||d t c„(i)|| 2 + — f \\d t r m (s)f ds + ^- I \\d t r m {s)-nf da ds 

UJ - Jo w - Jo 

< c(||a tC „(o)|| 2 + ||a t /|| 2 2(0iT;L2(O)) + \\d tg \\ 2 L2{w2{dn)) ) + / \\d tCn ( s ) \\ 2 d s , (3.9) 

Jo 

1 K 2 

with C = max(l, —5-, — — ). To bound ||<9tc n (0)|j, we use the first equation of (|3.3|) 

UJ_ UJ-Kl 

at t = with fi — dtc n , and the Cauchy-Schwarz inequality to obtain 

||ftc„(0)|| 2 < C(||/(0)|| 2 + ||V - r m (0)|| 2 ) < C(||/(0)|| 2 + llcoll^o))- 

Here we have used the fact that D~ 1 r m (0) = — Vc„(0) in 2?'(f2) given by the second 
equation of (|3.3|l . and hence in L 2 (Q) since Co £ H 2 (fl). From this inequality and 
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we have 

\\d t c n (t)\\ 2 + — [ \\d t r m (s)fds + ^- [ \\d t r m (s)-nf dn ds 

UJ - Jo w - Jo 

< C(ll c o|lH2 (a) + ||/||^i (0iT;i 2 (n)) + ||3l|ffi (0 .T;L2 (a n))) + / \\d t c n \\ 2 ds. (3.10) 

Jo 

It now follows from ()3.10ll and Gronwall's lemma that 

||9 t c„|| L ~(o,T;L 2 (o)) + \\9tr m \\ L 2^ T . L 2^ + \\d t r m • n|| L 2( 0T;i 2 (ao)) 

<C(II c o\\h 2 (U) + ll/lli/ 1 (0 : T;L2(n)) + llfl , ||ffi(0,T ; Z, 3 (an)) 

). (3.11) 



To obtain the estimate in the H(div 1 f2)-norm, we follow the same steps as for (J3Z 
to obtain 

|| V ' ?'m|lioc(o i T;L 2 (0)) - ^(\\ c o\\ 2 H 2 (fl) + ll/llff 1 (0,T;L 2 (n)) + llfl , lliyi(0,T;i 2 (af2)))' 

This along with the inequality (i) of Lemma 13.21 gives 

< ecu c o\\H 2 (n) + ll/llff 1 (0,T;L2(n)) + \\9\\H 1 (0,T-L 2 (dQ)) )• (3.12) 

The lemma now follows from (|3.11j) and (|3.12p . □ 

Thanks to Lemma 13.21 we can finish the proof of Theorem 13.11 using similar 
arguments as for the proof of Theorem 12.11 □ 

4. Space-time domain decomposition methods. In this section, we present 
two nonover lapping domain decomposition methods for solving problem (|2.3|) . For 
simplicity, we consider a decomposition of f2 into two non overlapping subdomains f^i 
and f^2 separated by an interface T: 

ii 1 nn 2 = Q>; r = afiinao 2 no, H^UiU^ur. 

Also for the sake of simplicity we have assumed throughout this section and the next 
that the boundary condition given on dQ is a homogeneous Dirichlet condition. How- 
ever, the analysis given below can be generalized to the case of multiple subdomains 
and more general boundary conditions (see Section [5]) . 

For i = 1,2, let denote the unit outward pointing vector field on <9fl;, and 
for any scalar, vector or tensor valued function tp defined on f2, let ip^ denote the 
restriction of ip to f2j. Using this notation, problem (|2.3p can be reformulated as 
an equivalent multidomain problem consisting of the following space-time subdomain 
problems 

w i 5 t c l + V-r l =/ inf2jX(0,T), 
Vc, + D^n = in n 4 x (0, T), _ 

Cl =0 on nan x (0,T), Iorz - i ' z ' 
Ci(0) = c in Oi, 

together with the transmission conditions on the space-time interface 

Cl = C \ n on T X (0, T) , (4.2) 
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Alternatively, and equivalently, one may impose the transmission conditions 

-ri-»i+ai 2 ci = -r 2 • Tlx + ai 2 c 2 r> „ m t\ (a i\ 

on T x (0,TJ , (4.3) 

-r 2 • n 2 + «2.i c 2 = -ri • n 2 + a 2 ,i c i 

where cti, 2 and a 2 .i are a pair of positive parameters. The first method that we 
consider is based on (|4.1[) together with the "natural" transmission conditions (|4.2[) 
while the second method is based on (|4.ip together with the Robin transmission 
conditions (|4.3[) . For the latter method the parameters ctij may be optimized to 
improve the convergence rate of the iterative scheme (see (T| ITO l \TT \ I12|). 

For both methods the multidomain problem is formulated through the use of in- 
terface operators as a problem posed on the space-time interface. For the first method 
the interface operators are time-dependent Steklov-Poincare (Dirichlet-to-Neumann) 
operators while for the second they are Robin-to-Robin operators. Associated with 
a Jacobi algorithm this latter method is known as the Optimized Schwarz Waveform 
Relaxation (OSWR) method. Rewriting the OSWR method as a space-time interface 
problem solved by a more general (Krylov) method was done in [H]; here we extend 
that work to a problem written in mixed form. 

4.1. Method 1: Using the time-dependent Steklov-Poincare operator. 

To introduce the interface problem for this method we introduce several operators, 
but first we define some notation: 

A = H 1 ^, T; hI q (T)), and, for i = 1, 2, M 4 = L 2 (Qi) and E 4 = ff(div, flij. 

We also define if*(fi,) = {v € H 1 ^), v = over dQ., n dQ}, fori = 1, 2. 
Next, let T>i, i = 1,2, be the solution operator that associates to the boundary, right- 
hand-side, and initial data (A, /, cq) the solution (cj,7"i) of the subdomain problem 

oJidta + V-n = f mfiiX(0,T), 

Vci+irVi =0 in fij x (0, T) , 

Cl = on dSk ndrix (0, T), (4.4) 

Cl = A on T x (0, T) , 

c,(0) = c inQi. 

An extension of Theorem 12.11 (to the case of non-homogeneous Dirichlet boundary 
conditions) guarantees that 

V, : Ax L 2 (0,T;L 2 (fiO) X #*(^) -> ^(0, T; Mi) x L 2 (0, T; E 4 ) 

(A,/, co) H' (ci,r<) = (c i (A,/,co),ri(A,/,c )) 

is a well defined operator. We also make use of the normal trace operator 

Ti : F 1 (0,T;M i ))xL 2 (0,T;E i ) — > L 2 (0, T; (ij| (r))') 
(ci.n) i y n -ni|rx(o,T) 

which is then used to define the following operators: 

5, : A — ► L 2 (0,T;(i/| (r))') 
A (-> (A, 0,0) 

and 

Xl : I 2 (0 1 T;I 2 (f! i ))xF, 1 (!l i ) -> L 2 (0,T;(i4,(r))') 
(/,c ) TO(0,/,c o ). 
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Now letting S = Si + S2 and \ = Xi + X2 we may rewrite problem (|4.ip , 2[) as the 
interface problem 

5A = x(/,co), onrx(0,T). (4.5) 
The weak formulation of this problem is then 

Find A € A such that: 

(4.6) 



(SX,r,) = / (x(/,co),»?>, \fi] e A, 
Jo 

where (•, •) denotes the duality pairing between Hq (T) and (iJ 2 (r))'. The operator 
S is the time-dependent Steklov-Poincare operator, and to investigate its properties 
we write the weak formulation of the interface problem (|4.4p for / = and cq = 0: 



For a.e. t £ (0,T), find a(t) <E M t and r,(t) G Si such that 
d_ 
dt 



— (uJiCi, fj,) Qi + (V • n, (j,) Qi = 0, V/U e Mi, (4.7) 



-(V + (Z^ Vi,t>) n< =-J\(v-ni), VweE,. 

For A 6 A and for i = 1,2, we will denote by (cs(A), (A)) the solution of (|4.7[) for the 
data function A. Then for 77, A g A and for almost every i S (0, T), we have 

(w 4 9 t Ci(A),Ci(ry))n ! + (V ■ rj(A), ^(17))^ =0, 

-(VTi^CifA))^ + (i) i - 1 r i (A),r i (77))n i = - / Hniv) • «*)• 



Now adding the first equation to the second equation in which the roles of A and 77 
are reversed, integrating over time and summing on i, we obtain 

V/ ((^a tCi (A), c^))^ +(D?tM, r,(A)) 0l ) =-E/ fviriW-ni). 
l=1 Jo i=1 Jo Jr 



Thus we see that 
r T 2 „t 



/ (SX,r,) =-J2 [ f (r<(A) = V / ((^^(A), c^))^ + (DrV^A),^))^) , 
./o 4=1 ./o./r t=1 Jo 

from which we conclude that S is a positive definite but non-symmetric, space-time 
interface operator. Thus the existence and uniqueness of the solution of the space-time 
interface problem (|4.6[) does not follow in a standard way, and we have not pursued 
this question here. 

Nonetheless, we solve a discretized version of problem (|4.5j) iteratively by using 
a Krylov method (e.g. GMRES). Once the discrete approximation to A is obtained, 
we can construct the multi-domain solution of the discretized problem. Following 
the work in [22l [25] for elliptic problems with strong heterogeneities, we apply a 
Neumann-Neumann type preconditioner enhanced with averaging weights: 

(a 1 5r 1 +a 2 5 2 - 1 )5A = x, (4.8) 

where at : T x (0,T) — > [0,1] is such that <j\ + o"2 = 1, and S^ 1 , the Neumann-to- 
Dirichlet operator, is the (pseudo)-inverse of <Sj, for i = 1,2. 
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4.2. Method 2: Using Optimized Schwarz Waveform Relaxation (OSWR). 

The function spaces that are needed to give the interface formulation of method 2 are 

S:=iT 1 (0,T;L 2 (r)) ! and, for i = 1, 2, Af, = L 2 (fi,) and Si = H(div, Qj). 

To define the Robin-to- Robin operator we first define for i = 1,2, the following solution 
operator TZi which depends on the parameter Oiij; j = 3 — i : 

Hi-. SxL^O.T^^Dx^f!!,) — > 3 x H^O^M,) x L 2 (0,T;£,) 

(£,/,c ) i— > Ci,ri) = (^,Ci(C,/,co),ri(^,/, co)) 

where (cj,rj) = (c,(£, /, Co),t"i(£, /, Co)) is the solution to the problem 

Wi9 t c l + V • r l = / in f2i x (0, T) , 
Vci + DrV* =0 in f2j x (0, T) , 

c t =0 on dCli Ddflx (0, T), (4.9) 
-Ti ■ rii + otijCi =i on T x (0, T) , 
c,(0) = c infij. 

(As stated earlier the parameters ajj will be chosen is such a way as to optimize 
the convergence of the algorithm). The existence and uniqueness of the solution of 
problem (|4.9|) is guaranteed by Theorem 13. II 

Next, to impose the interface conditions (|4.3[) we will need the following interface 
operators defined for i = 1,2, and j — 3 — i: 

Bi'. (e x H\0,T;Mj) x L 2 (0,T;£ i )') nlm(^) — > 3 

(£.c,-,rj) >-> (-rj •n i + -^-(^ + r i • n.,-))|rx(o,T) 

Remark 4.1. To see t/iat Im(Bi) C H (instead of simply L 2 (0,T; L 2 (r)), we noie 
tfiaf d3T2]) implies t/iat i) -1 r(i) = -Vc(i) in X>'(fi) for a.e. t £ (0,T). Since r(t) is 
in H(div,Cl), we have c(t) £ H (CI), for a.e. t £ (0,T). Consequently, c,(t) is in 
-ff^O, T; if (Oj)). XTiis along with the fact that implies that • «i|rx(o,T) € 2. 



Now, defining 



S H : 3 x 3 — > ^ x 3 

/a\ /a -si^ 2 (6,o,o) 
UJ ^ 16-^1 (a, o,o) 



and 



XR : L 2 (0,r;i 2 (a)) x Hl{fli) 
(f,co) 

we can write the interface problem as 



£ift 2 (0,/, c ) 
B 2 Ki{0,f, c ) 



5 *(|) =Xft(/,Co) onrx(0,T). (4.10) 
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We then write (|4.10p in weak form as 
Find ((i,( 2 )6HxH such that 

.>-(!) ■(D-n^'-Cc;)-^ 6 ' 63 ^ <4,11) 



In order to study the interface operator Sr, we proceed as in the Section [4.1l bv giving 
the weak formulation of the relevant subdomain problems (here ()4.9|) for i = 1, 2 and 

j = 3 — 1) for / = and Co = 0: 

For a.e. t G (0,T), find Cj(i) € Mi and r*(i) G such that, V> G M t and Vv G E;, 

— {u} i c i ,n)Q i + (V • r 4 , /Lt)n s = 0, 

-(V -WjCi)^ + (•Dr lr t,«)n< + / (*•» -ni)(v -n,) = - / — £(«-nj). 

Jr a i,j Jr a i,j 

(4.12) 

Now for any ( G S letting a{Q G ff^O, T; M») and n(() G L 2 (0, T; Si) be such that 
7£,(C, 0,0) = (C)Cj(C)) r i(C))i we have for any pair of elements £ and £ in S and for 
a.e. t G (0, T) that 

(uidtCt(t)MQ)n< + (V- CiCOJn, = 0, 
- (V • r< (C) , Q (Ok + (Dr 1 * 1 * (0 , r, (OK + / — fa (0 • n*) fa (0 ■ * ) 

— £fa(C)-»i)- 

r o»,3 

Next we add the first of these two equations to the second in which the roles of £ and 
£ have been interchanged to obtain 

fadtCiiQMOot + (D^niOMOh, + f — (r,(0 -nOfate) -m) 

jr - - (4.13) 



/ — C(r<(0-«< 



'r 

and this holds for any pair of elements £ and £ in E. Now we consider the case in 
which the parameters on t j, i = 1,2, j = 3 — i, are constant and apply (|4.13p with 
£ = £j and C = Ci; to obtain 

Jo L Sr (&) ' (^) = " 57 G ' )Ci + K2 + a2 ' 1 ){^ 9tC '^)' c '^)K 

+(£>r 1 ^fe),rt(0))n i + / — fa(&) - Hi)}} 

As for method 1, we obtain a non-symmetric, space-time interface operator, but here it 
is also not positive definite. We solve the discretized problem iteratively using Jacobi 
iterations or GMRES. The former choice is equivalent to the OSWR algorithm, and 
in the next subsection we show that this mixed form of the algorithm converges. 

4.2.1. The OSWR algorithm. We consider the general case in which £1 is 
decomposed into I non-overlapping subdomains f2j. We denote by Fj j the interface 



21 



between two neighboring subdomains Oj and Oj, Tjj = dili n 9fij n fi. Let Mi be 
the set of indices of the neighbors of the subdomain Oj, i = 1, . . . , J. The OSWR 
method may be written as follows: at the fc" 1 iteration, we solve in each subdomain 
the problem 

d t c k + V-r k = /, mJ],x(0,r), 

Vcf + DrVf = 0, in fii x (0, T) , (4.14) 

-r\ -rii+aijc^ = -r^ 1 ■ n t + a iijC k r l , onT^ x (0,T),Vj e M, 



where, for i = 1, • • • , I, j £ Mi, otij > is a Robin parameter. The initial value is 



that of Co in each subdomain. Moreover, (gi.j) := — r° • rii + otijc® is an initial guess 



on IYj, for i = 1,- • ■ Mi, in order to start the first iterate. 

Theorem 4.2. Lef £> € W^ 1 ' 00 ^), / e /^(O^L 2 ^)) and c £ H 2 (fl)nH^(n) 
and let ai.j € L co (dfli) be such that ajj > ao > /or i = 1, ••■ j e A^. 
Algorithm (|4.14p . initialized by (gij) in H 1 (0, T; L 2 (Tjj)) , i = 1, • • ■ J € A/i, 
defines a sequence of iterates 

(^rf) eW^O^L 2 ^)) x (L 2 (0,T;H(div,^))ni/ 1 (0,T;L 2 (n i ))), 

/or i = 1, •• • , 7 ; £/ia£ converges to the weak solution (c,r) of problem (|2.3p . 

Proof. The sequence (c^,r^)fe is well-defined according to Theorem 13.11 and Re- 
mark [4J] Now, to prove the convergence of algorithm (|4.14|1 . as the equations are 
linear, we can take / = and cq ~ and show that the sequence [c k ,r k ) k of iterates 
converges to zero in suitable norms. 

To begin, we write the variational formulation of (|4.14p (with / = 0): 



For a.e. t e (0, T). find c k (t) £ M t and r k (t) G Sj such that 

-(wcf./ii)^ + (V-r, fc ,^)n s =0, V[m £ M h 

-{V ■v i ,c k )n, + {Df 1 r k ,v i )n, = E / c i(~ v i • 0> V»« 6 ^. 



(4.15) 



Choosing /i.; = and U; = in (|4.15[) , then adding the two resulting equations and 
replacing the boundary term by using the equation 



■ n + a Z iC k ) - (-r k ■ n, - a.-.-c* 



2 (aij + ccj,i) c k (~r k ■ m) + (al 3 - a\ ,) (c k ) 



we obtain 



We then integrate over (0,i) for a.e. i € (0,T] and apply the Robin boundary con- 
ditions. By using the properties of uj and D and recalling that the Robin coefficients 
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ctij belong to L°° i S 1, • • • , I, j S Mi, we obtain, for some constant C, 
W _||c?(t)||^+M_ / t ||rf( S )|| 2 ds+ E T / T ( r 1 ; n, ,,,.,,>.) 

Jo jT iti a h] + a 3^ JO 

Now we sum over all subdomains and define for k > 1 and for a.e. t € (0, T] 
E k (t) = J2 (w_||cf (t) |& + 2<5_ £ ||r? (a) \\ 2 ds 



i=l 
7 



Then we have, for all k > 

7 /•« 



E k (t) + B k (t) < B k ^ (t) + C^ [ \\c k (s) f Qi ds. 



Now sum over the iterates for any given K > 

K K I t 



j2E k (t)<B°(t) + cj2J2 1 ii c "( s )ii^ ( 4 - 16 ) 

fc=l k=l i=l ^° 



where 



for gij the initial guess on Tij. From the definition of E , since J_ > 0, we have 

E I>-ii c * (*) n* < 5° (*) + ^EE T ii c " w n**- 

fe=l i=l fe=l i=l ^° 

Thus, by applying Gronwall's lemma, we obtain for any K > and a.e. t £ (0,T) 

K I 



CT 



fe=l i=l 

This along with (|4.16p implies 

K I 



EE 2 ^ \\rHs)\\ 2 d S <(l + ^—e-)B a (T), VK > 0. (4.18) 



k=l i=l 



The inequalities (|4.17p , (|4.18p imply that the sequence c\ tends to in L°° (0, T; L 2 (fi*)) 
and converges to in L 2 (0, T; X 2 (f2j)) for each i £ 1, • • ■ ,/asfc— ^ oo. 
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To show convergence in higher norms, we differentiate the first and the second equa- 
tions of f|4. 1 5[) with respect to t, then take fii — dtc k and Vi = dtr k and add the 
resulting equations together, we see that (after bounding the left hand side using the 
assumptions on oj and D) 



ui_ d 

~T~dr 



d t c k \\ 2 Qi + 6-WdtrUk < E / d ^(-9tr1 ■ 



We proceed as in the previous argument with the use of Robin boundary conditions 
after differentiating with respect to t 

-Serf ■ n, + <>,.Jh< } ; = -dt^f 1 ■ n t + a^dtc) : . on T id x (0, T) , Vj € M- 

We then obtain, for all k > 

E k (t) + B k (t) < B k ~ 1 (t) +C V / ||9 t c* (s) ||^d«. 

where 



^ (*) = E {^-W 9 ^ (*) I& + 2*- f Q \\dtr k (a) \\ 2 ds 

B k (*) = EE /*/ ^4^- (-^ • * + "<-,'Kr . 

Now, as before, we sum over the iterates for any K > and apply Gronwall's lemma 
to obtain for any K > and a.e. t £ (0, T) 



EEii^wii^< e ^^, with b° (*) = E E f / — 1 — (^) 2 

(4.19) 

This along with (|4~T7| shows that the sequence c k converges to in W l oo (Q, T; L 2 (f2i) 
as k — > oo, for i = 1, 

Now we choose /Zj = V • in the first equation of (|4.15|) to obtain for a.e. t € (0, T) 
||V-r*|| 2 = V-r*) <||S t cf||||V-^||. 

or 

llV-rfU <||9 t ^|| We(0,T). 



Hence, by (|4.19j) we have 

||V • r k \\ Lac{0tT . L 2 {nt)) -> as fc -> oo. (4.20) 

This shows that the sequence r k converges to in L 2 (0, T; i?(div, fli)). Moreover, it 
follows from the definition of E k and (|4.19jl that 

K 1 ft (~lrp CT 

EE 25 - / \\dtr k (s)\\ 2 ds<(l + —e~)B Q (T), MK > 0. 

Jfc=l i=l ^° W ~ 

So that the sequence d t r k also converges to in L 2 (0, T; L 2 (fli)). □ 
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5. Nonconforming time discretizations and projections in time. One of 

the main advantages of Method 1 or Method 2 is that these methods are global in 
time and thus enable the use of independent time discretizations in the subdomains. 
At the space-time interface, data is transferred from one space-time subdomain to a 
neighboring subdomain by using a suitable projection. 

We consider semi-discrete problems in time with nonconforming time grids. Let 71 
and 72 be two possibly different partitions of the time interval (0, T) into sub-intervals 
(see Figure ISTTj) . We denote by the time interval (t t m _ 1 ,t t m ] and by At l m := (tf m — 
t l m -i) f° r Tin = 1, ... , Mi and i = 1,2, where for simplicity of exposition we have again 
supposed that we have only two subdomains. We use the lowest order discontinuous 
Galerkin method [3l [TTl [33] , which is a modified backward Euler method. The same 
idea can be generalized to higher order methods. We denote by Po(77, W) the space of 



T 



Atl 



At 2 



Qi n 2 x 

Figure 5.1: Nonconforming time grids in the subdomains. 

piecewise constant functions in time on grid % with values in W, where W = (T) 
for Method 1 and W = L 2 (T) for Method 2: 

P (%, W) = {<f> : (0, T) -> W, <f> is constant on 4, Mm = 1, . . . , Mi} . 

In order to exchange data on the space-time interface between different time grids, we 
define the following L 2 projection IF,, from P (Ti,W) onto P (Tj,W) (see [HlfTT]) : 
for <p € Po{Ti, W), Hji4>\jj is the average value of <fr on J ? 3 n , for m = 1, . . . , Mj: 

Mi 

\ Jfn \^~[ JJLnJ; 

We use the algorithm described in [T3] for effectively performing this projection. With 
these tools, we are now able to weakly enforce the transmission conditions over the 
time intervals. 

We still denote by (ci,r;), for i = 1,2, the solution of the problem semi-discrete in 
time corresponding to problem (|4.7I) or (|4.12D . 

5.1. For Method 1. As there is only one unknown A on the interface, we need 
to choose A piecewise constant in time on one grid, either 7i or 72- For instance, 
let A G P (72, -ff^(r)) and take c 2 = n 22 (A) = Id(A). The weak continuity of the 
concentration in time across the interface is fulfilled by letting 

ci=n 12 (A) g P {Ti,Hi(T)). 

The semi-discrete (nonconforming in time) counterpart of the flux continuity in the 
second equation of (|4.2[) is weakly enforced by integrating it over each time interval 
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J m of grid 72 : Vm = 1, M 2 , 

n 21 (n (n 12 (A), /, c ) • m) + n 22 (r 2 (n 22 (A), /, c ) • n 2 ) ] dt = o. (5.1) 



Remark. Obviously one can choose A to be constant in time on yet another grid 
(neither 71 nor T2), and this can be useful in some applications (e.g. flow in porous 
media with fractures). 

5.2. For Method 2. In Method 2, there are two interface unknowns represent- 
ing the Robin terms from each subdomain. Thus we let £3 G Po{7~i, L 2 (T)), for i = 1,2. 
The semi-discrete in time counterpart of (|4.3|) is weakly enforced as follows: 



r jji 



ii - IIi 2 (-r 2 (£ 2! /,co) Tlx + ai, 2 c 2 (£ 2 ,/,c ))J dt = 0, Vm = 1,- • • ,Mi, 

-n 2 i(-ri(^i,/,c ) • n 2 + a 2 ,ici(^i, /, c )) +6^ * = 0, Vm = 1, • • • ,M 2 , 

(5.2) 



/r J 

where (ci(&, /, c ), /, c )) , i = 1,2 is the solution to (|4.12|) . 



Remark. For conforming time grids, the two schemes defined by applying GMRES for 
the two interface problems (|5.1|) , (|5.2p respectively converge to the same monodomain 
solution. In the nonconforming case, due to different projections, the two schemes 
become different and in the next section, we will study and compare the errors in time 
for the two approaches. 

6. Numerical results. In this section, we carry out numerical experiments in 
2D to illustrate the performance of the two methods presented above. We consider 
D = dl isotropic and constant on each subdomain, where I is the 2D identity matrix. 
Consequently, we may denote by di, the diffusion coefficient in the subdomains. For 
the spatial discretization, we use mixed finite elements with the lowest order Raviart- 
Thomas spaces on rectangles [51 130]. 

In the first test problem (see Section l6.1j) . we consider the two subdomain case 
with discontinuous coefficients. We vary the jumps in the diffusion coefficients and 
we see how it affects the convergence speed. We also analyze the behavior of the error 
versus the time steps in the nonconforming case. In the second test problem (see 
Section l6T2|) . suggested by ANDRA as a first step towards repository simulations, we 
consider several subdomains. We observe how both methods handle this application 
with the strong heterogeneity and long time computations. 

6.1. A two subdomain case. The computational domain f2 is the unit square, 
and the final time is T = 1. We split into two nonover lapping subdomains 
Oi = (0,0.5) x (0, 1) and fi 2 = (0.5, 1) x (0, 1) as depicted in Figure IO The initial 
condition is cq — cxp ((x - 0.55) 2 + 0.5(y - 0.5) 2 ) and the right-hand side is / = 0. 
The porosity is uii = u> 2 = 1, the diffusion coefficients are d\ and c? 2 in Q 1 and f2 2 
respectively (d\ 7^ c? 2 ). We fix d 2 = 0.2 and vary d\ as shown in Table RTTI We let 
J) denote the diffusion ratio <i 2 /di. For the spatial discretization, we use a uniform 
rectangular mesh with size Axi = Ax 2 = 1/200. For the time discretization, we use 
nonconforming time grids with Ati and At 2 , given in Table [6T| adapted to different 
diffusion ratios. We first analyze the convergence behavior of each method. We solve 
a problem with cq = and / = (thus c = and r = 0). We start with a random 
initial guess on the space-time interface. We remark that one iteration of Method 1 
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Figure 6.1: Domain decomposition and boundary conditions. 
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d 2 


1/At 2 


10 


0.02 


150 


0.2 


200 


100 


0.002 


50 


0.2 


200 


1000 


0.0002 


20 


0.2 


200 



Table 6.1: Diffusion coefficients and corresponding nonconforming time steps. 



with the preconditioner costs twice as much as one iteration of Method 2 (in terms 
of number of subdomain solves). Thus to compare the two approaches, we plot the 
error (in logarithmic scale) in the L 2 (0, T; L 2 (i7))-norm of the concentration c and the 
vector field r, versus the number of subdomain solves (instead of versus the number 
of iterations) . We stop the iteration when the errors (both in c and r) are less than 
10~ 6 . In Figure HT2"1 the convergence of the two methods (with GMRES) for different 
diffusion ratios is shown. We see that both methods work well. Method 1 (Schur) 




D = 10 D = 100 D = 1000 

Figure 6.2: Convergence curves for different diffusion ratios: errors in c for Method 1 
(red) and Method 2 (blue); errors in r for Method 1 (magenta) and 
Method 2 (green). 



converges faster than Method 2 (Schwarz) for small diffusion ratios 2). However, when 
£> is increased, they are comparable. We also observe that the errors in c and r are 
nearly the same for Method 2 while the error in r is greater than the error in c for 
Method 1. Both methods handle the heterogeneities efficiently. To obtain such a good 
performance, we have used the following formula for calculating the weights in (14.81) 
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(see [23]) 

Consider now the case with D = 10. For Method 2, we vary Robin parameters ai,2 
and CH2.1 and plot the logarithmic scale of the residual after 20 Jacobi iterations in 
Figure [6731 We see that the optimized Robin parameters (the red star), which are 
calculated by numerically minimizing the convergence factor [TJ [2j 111) , are located 
close to those giving the smallest residual after the same number of iterations. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 6.3: Level curves for the residual (in logarithmic scale) after 20 Jacobi 
iterations for various values of the parameters ot\.i and 0.2,1- The red star shows the 
optimized parameters computed by numerically minimizing the continuous 

convergence factor. 

Next, we analyze the accuracy in time for different diffusion ratios and corre- 
sponding choices of nonconforming time steps. Toward this end, we consider four 
initial time grids (for At c and Atf given) 

• Time grid 1 (fine- fine): conforming with Ati = At2 = Atf. 

• Time grid 2 (coarse-fine): nonconforming with Ati = At c and Ai2 = Atf. 

• Time grid 3 (fine-coarse): nonconforming with Ati = Atf and Ai2 = At c . 

• Time grid 4 (coarse-coarse): conforming with Ati = At2 = At c . 

The time steps are then refined several times by a factor of 2. In space, we fix a 
conforming rectangular mesh and we compute a reference solution by solving prob- 
lem (|2.4I) directly on a very fine time grid, with At = Atf/2 e . The converged mul- 
tidomain solution is such that the relative residual is smaller than 10 . We show 
in Figures IBT41 and |6~51 the errors in the L 2 (0, T; L 2 (f2))-norms of the concentration c 
and the vector field r versus the time step At — max(Ai c , Atf) for different diffusion 
ratios. We only give the results for Method 1 because the curves for Method 2 look 
exactly the same. For D = 10, we take At c = 1/94 and Atf = 1/128; for D = 100, 
we take At c = 1/40 and Atf = 1/160 (for D = 1000, the same results hold for 
At c = 1/16 and Atf = 1/160 but we don't present it here). We first observe that 
first order convergence is preserved in the nonconforming case. Moreover, the error 
obtained in the nonconforming case (Time grid 2, in blue) is nearly the same as in 
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the finer conforming case (Time grid 1, in red). This means that nonconforming time 
grids preserve the solution's accuracy in time and one must refine the time step where 
the solution varies most (i.e. where the diffusion coefficient is larger). 




Figure 6.4: Errors in c (left) and r (right) in logarithmic scales between the 
reference and the multidomain solutions versus the time step for 2D = 10. 




At At 

Figure 6.5: Errors in c (left) and r (right) in logarithmic scales between the 
reference and the multidomain solutions versus the time step for 2) = 100. 

6.2. A porous medium test case. In this subsection, we consider a simplified 
version of a problem simulating contaminant transport in and around a nuclear waste 
repository site. The test case is described in Figure [6~B"1 where the repository is shown 
in red and the clay layer in yellow. The domain is a 3950m by 140m rectangle and 
the repository is a centrally located 2950m by 10m rectangle. The initial condition is 
Co = 0, the source term is / = in the clay layer and 

, f lO^s- 1 if t < 10 5 years, . . 

/ = < _ , .-k m the repository. (6.1 

J \ if t > 10 5 years, y J 

We impose homogeneous Dirichlet conditions on top and bottom, and homogeneous 
Neumann conditions on the left and right hand sides. We decompose H into 9 sub- 
domains as depicted in Figure 16.71 with fis representing the repository. The porosity 
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is CJ5 = 0.2 and a;, = 0.05, i 7^ 5. The diffusion coefficients are d§ = 2f0- 9 m 2 /s 
and di = 5 10 -12 m 2 /s, i 7^ 5. So the diffusion ratio is ID = 400. For the spatial 



140 m 




3950 m 

Figure 6.6: Geometry of the domain. 
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Figure 6.7: The decomposition into 9 subdomains (blow up in the y-direction) . 

discretization, we use a non-uniform but conforming rectangular mesh with a finer 
discretization in the repository (a uniform mesh with 600 points in the x direction and 
30 points in the y direction) and a coarser discretization in the clay layer (the mesh 
size progressively increases with distance from the repository by a factor of 1.05). For 
the time discretization, we use nonconforming time grids with At^ = 2000 years and 
Ati = 10, 000 years, i 7^ 5. For this application, we are interested in the long-term be- 
haviour of the repository, say over one million years. Thus, we test the performance of 
the two methods for a "short" time interval (T = 200, 000 years) and for a longer time 
interval (T = 1,000,000 years). The same time steps, Ati, are used for both cases. 
As in the first test problem, we analyze the convergence results by solving a problem 
with / = 0. For Method 2, as we have a small, thin object embedded in a large area, 
it has been shown in |16[ 120) that it is important to derive an adapted optimization 
for Robin parameters. Thus, we consider two different optimization techniques: the 
classical one (Opt. 1) as used in the first test problem, and an adapted version (Opt. 2) 
[TBI |2"0] where we take into account the dimension of the subdomains. 

In Figure [6751 we compare the errors in the concentration c (on the left) and in the 
vector field r (on the right) both over a shorter time interval (on top) and over a longer 
time interval (on bottom) where GMRES is used in all cases as the iterative solver: 
Method 1 (red), Method 2 with Opt. 1 (blue) and Method 2 with Opt. 2 (green). They 
are comparable and perform well in the case of multiple subdomains. We also note 
that the longer the time interval, the larger the number of subdomain solves needed 
to converge to a given tolerance (here 10~ 6 ). Thus, the use of time windows (see 
[31 [TB]) could considerably improve the performance of all the algorithms, especially 
with an adapted choice of the initial guess on the interface based on the solution on 
the previous time window. In Figure 16.91 we plot the errors in the concentration c 
over different time intervals for Method 2 with Jacobi iteration: with Opt. 1 (blue) 
and Opt. 2 (green) (the errors in the vector field r behave similarly). We observe that 
Opt. 2 efficiently handles the long time computation case while Opt. 1 doesn't. 

Next we consider the case with / 7^ as defined in (|6.1|) and over the long time 
interval, T = 1, 000, 000 years. The discretizations in space and in time (nonconform- 
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Figure 6.8: Convergence curves for different time intervals with GMRES: error in c 
(on the left) and error in r (on the right), for short time T = 200, 000 years (on top) 
and for long time T = 1, 000, 000 years (on bottom). 

ing) are the same as above. We verify the performance of Method 1 and Method 2 
(with Opt. 2) using GMRES and zero initial guess on the space-time interfaces. The 
tolerance of the iteration is 10 -6 . In Figure IS. 101 the evolution of the solution at 
different times is depicted (both methods give similar results). As time goes on and 
under the effect of diffusion, the contaminant slowly migrates from the repository to 
the surrounding area. Moreover, its concentration c increases until injection stops (i.e. 
after 100,000 years) and then decreases. In Figure IB. Ill the relative residuals for each 
method versus the number of subdomain solves are shown, as the monodomain solu- 
tion with nonconforming grids is unknown. Both methods work well and we observe 
that Method 1 converges linearly while Method 2 initially converges extremely rapidly, 
the convergence becoming linear after the first few iterations. 

7. Conclusion. We have given mixed formulations for two different interface 
problems for the diffusion equation, one using the time-dependent Steklov-Poincare 
operator and the other using OSWR with Robin transmission conditions on the space- 
time interfaces between subdomains. The subdomain problem with Robin boundary 
conditions is proved to be well-posed. A convergence proof of the OSWR algorithm 
in mixed form is also given. Nonconforming time grids are considered and a suitable 
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Figure 6.9: Convergence curves for different time intervals using Jacobi iteration: for 
short time T = 200, 000 years (on the left) and for long time T = 1, 000, 000 years 

(on the right). 




Figure 6.10: Snapshots of the multi-domain solution after 20,000 years (top left), 
100 000 years (top right), 200 000 years (bottom left), and 1,000,000 years (bottom 
right), with a blow up in the y-direction. 

projection in time is employed to exchange information between subdomains on the 
space-time interface. Numerical results for 2D problems using mixed finite elements 
(with the lowest order Raviart-Thomas spaces on rectangles) for discretization in 
space and the lowest order discontinuous Galerkin method for discretization in time 
are presented. We have analyzed numerically the performance of the two methods for 
two test cases, one academic with two subdomains and one more realistic with sev- 
eral subdomains. We have observed that both methods handle well the heterogeneity 
and nonconforming time grids, both efficiently preserving the solution's accuracy in 
time. The two methods are also well-adapted for the simulation of diffusive contam- 
inant transport in and around a repository with a special geometry and long time 
computations. In particular, for Method 2 we have shown that an adapted optimiza- 
tion technique to compute the optimized parameters is necessary if Jacobi iteration is 
used. We have pointed out the possible advantage for efficiency of using time windows 



32 




10 15 20 25 30 35 '» „ g ,„ , 5 2 „ 25 

Number of Subdomain Solves Number of Subdomain Solves 



Figure 6.11: The relative residuals in logarithmic scales using GMRES for Method 1 
(on the left) and Method 2 (with Opt. 2) (on the right). 

for problems with long time interval. Work underway addresses the coupling between 
advection and diffusion using operator splitting as well as nonmatching grids in space. 
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